<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>shsein.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="SHSEIN.1"></a><a href="shsein.f.html#SHSEIN.1">SHSEIN</a>( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI,
     $                   VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL,
     $                   IFAILR, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          EIGSRC, INITV, SIDE
      INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      LOGICAL            SELECT( * )
      INTEGER            IFAILL( * ), IFAILR( * )
      REAL               H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
     $                   WI( * ), WORK( * ), WR( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="SHSEIN.23"></a><a href="shsein.f.html#SHSEIN.1">SHSEIN</a> uses inverse iteration to find specified right and/or left
</span><span class="comment">*</span><span class="comment">  eigenvectors of a real upper Hessenberg matrix H.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The right eigenvector x and the left eigenvector y of the matrix H
</span><span class="comment">*</span><span class="comment">  corresponding to an eigenvalue w are defined by:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">               H * x = w * x,     y**h * H = w * y**h
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where y**h denotes the conjugate transpose of the vector y.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  SIDE    (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'R': compute right eigenvectors only;
</span><span class="comment">*</span><span class="comment">          = 'L': compute left eigenvectors only;
</span><span class="comment">*</span><span class="comment">          = 'B': compute both right and left eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  EIGSRC  (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          Specifies the source of eigenvalues supplied in (WR,WI):
</span><span class="comment">*</span><span class="comment">          = 'Q': the eigenvalues were found using <a name="SHSEQR.43"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>; thus, if
</span><span class="comment">*</span><span class="comment">                 H has zero subdiagonal elements, and so is
</span><span class="comment">*</span><span class="comment">                 block-triangular, then the j-th eigenvalue can be
</span><span class="comment">*</span><span class="comment">                 assumed to be an eigenvalue of the block containing
</span><span class="comment">*</span><span class="comment">                 the j-th row/column.  This property allows <a name="SHSEIN.47"></a><a href="shsein.f.html#SHSEIN.1">SHSEIN</a> to
</span><span class="comment">*</span><span class="comment">                 perform inverse iteration on just one diagonal block.
</span><span class="comment">*</span><span class="comment">          = 'N': no assumptions are made on the correspondence
</span><span class="comment">*</span><span class="comment">                 between eigenvalues and diagonal blocks.  In this
</span><span class="comment">*</span><span class="comment">                 case, <a name="SHSEIN.51"></a><a href="shsein.f.html#SHSEIN.1">SHSEIN</a> must always perform inverse iteration
</span><span class="comment">*</span><span class="comment">                 using the whole matrix H.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INITV   (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'N': no initial vectors are supplied;
</span><span class="comment">*</span><span class="comment">          = 'U': user-supplied initial vectors are stored in the arrays
</span><span class="comment">*</span><span class="comment">                 VL and/or VR.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  SELECT  (input/output) LOGICAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          Specifies the eigenvectors to be computed. To select the
</span><span class="comment">*</span><span class="comment">          real eigenvector corresponding to a real eigenvalue WR(j),
</span><span class="comment">*</span><span class="comment">          SELECT(j) must be set to .TRUE.. To select the complex
</span><span class="comment">*</span><span class="comment">          eigenvector corresponding to a complex eigenvalue
</span><span class="comment">*</span><span class="comment">          (WR(j),WI(j)), with complex conjugate (WR(j+1),WI(j+1)),
</span><span class="comment">*</span><span class="comment">          either SELECT(j) or SELECT(j+1) or both must be set to
</span><span class="comment">*</span><span class="comment">          .TRUE.; then on exit SELECT(j) is .TRUE. and SELECT(j+1) is
</span><span class="comment">*</span><span class="comment">          .FALSE..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix H.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  H       (input) REAL array, dimension (LDH,N)
</span><span class="comment">*</span><span class="comment">          The upper Hessenberg matrix H.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDH     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array H.  LDH &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WR      (input/output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">  WI      (input) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, the real and imaginary parts of the eigenvalues of
</span><span class="comment">*</span><span class="comment">          H; a complex conjugate pair of eigenvalues must be stored in
</span><span class="comment">*</span><span class="comment">          consecutive elements of WR and WI.
</span><span class="comment">*</span><span class="comment">          On exit, WR may have been altered since close eigenvalues
</span><span class="comment">*</span><span class="comment">          are perturbed slightly in searching for independent
</span><span class="comment">*</span><span class="comment">          eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  VL      (input/output) REAL array, dimension (LDVL,MM)
</span><span class="comment">*</span><span class="comment">          On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
</span><span class="comment">*</span><span class="comment">          contain starting vectors for the inverse iteration for the
</span><span class="comment">*</span><span class="comment">          left eigenvectors; the starting vector for each eigenvector
</span><span class="comment">*</span><span class="comment">          must be in the same column(s) in which the eigenvector will
</span><span class="comment">*</span><span class="comment">          be stored.
</span><span class="comment">*</span><span class="comment">          On exit, if SIDE = 'L' or 'B', the left eigenvectors
</span><span class="comment">*</span><span class="comment">          specified by SELECT will be stored consecutively in the
</span><span class="comment">*</span><span class="comment">          columns of VL, in the same order as their eigenvalues. A
</span><span class="comment">*</span><span class="comment">          complex eigenvector corresponding to a complex eigenvalue is
</span><span class="comment">*</span><span class="comment">          stored in two consecutive columns, the first holding the real
</span><span class="comment">*</span><span class="comment">          part and the second the imaginary part.
</span><span class="comment">*</span><span class="comment">          If SIDE = 'R', VL is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDVL    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array VL.
</span><span class="comment">*</span><span class="comment">          LDVL &gt;= max(1,N) if SIDE = 'L' or 'B'; LDVL &gt;= 1 otherwise.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  VR      (input/output) REAL array, dimension (LDVR,MM)
</span><span class="comment">*</span><span class="comment">          On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
</span><span class="comment">*</span><span class="comment">          contain starting vectors for the inverse iteration for the
</span><span class="comment">*</span><span class="comment">          right eigenvectors; the starting vector for each eigenvector
</span><span class="comment">*</span><span class="comment">          must be in the same column(s) in which the eigenvector will
</span><span class="comment">*</span><span class="comment">          be stored.
</span><span class="comment">*</span><span class="comment">          On exit, if SIDE = 'R' or 'B', the right eigenvectors
</span><span class="comment">*</span><span class="comment">          specified by SELECT will be stored consecutively in the
</span><span class="comment">*</span><span class="comment">          columns of VR, in the same order as their eigenvalues. A
</span><span class="comment">*</span><span class="comment">          complex eigenvector corresponding to a complex eigenvalue is
</span><span class="comment">*</span><span class="comment">          stored in two consecutive columns, the first holding the real
</span><span class="comment">*</span><span class="comment">          part and the second the imaginary part.
</span><span class="comment">*</span><span class="comment">          If SIDE = 'L', VR is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDVR    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array VR.
</span><span class="comment">*</span><span class="comment">          LDVR &gt;= max(1,N) if SIDE = 'R' or 'B'; LDVR &gt;= 1 otherwise.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  MM      (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of columns in the arrays VL and/or VR. MM &gt;= M.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  M       (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of columns in the arrays VL and/or VR required to
</span><span class="comment">*</span><span class="comment">          store the eigenvectors; each selected real eigenvector
</span><span class="comment">*</span><span class="comment">          occupies one column and each selected complex eigenvector
</span><span class="comment">*</span><span class="comment">          occupies two columns.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) REAL array, dimension ((N+2)*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IFAILL  (output) INTEGER array, dimension (MM)
</span><span class="comment">*</span><span class="comment">          If SIDE = 'L' or 'B', IFAILL(i) = j &gt; 0 if the left
</span><span class="comment">*</span><span class="comment">          eigenvector in the i-th column of VL (corresponding to the
</span><span class="comment">*</span><span class="comment">          eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
</span><span class="comment">*</span><span class="comment">          eigenvector converged satisfactorily. If the i-th and (i+1)th
</span><span class="comment">*</span><span class="comment">          columns of VL hold a complex eigenvector, then IFAILL(i) and
</span><span class="comment">*</span><span class="comment">          IFAILL(i+1) are set to the same value.
</span><span class="comment">*</span><span class="comment">          If SIDE = 'R', IFAILL is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IFAILR  (output) INTEGER array, dimension (MM)
</span><span class="comment">*</span><span class="comment">          If SIDE = 'R' or 'B', IFAILR(i) = j &gt; 0 if the right
</span><span class="comment">*</span><span class="comment">          eigenvector in the i-th column of VR (corresponding to the
</span><span class="comment">*</span><span class="comment">          eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
</span><span class="comment">*</span><span class="comment">          eigenvector converged satisfactorily. If the i-th and (i+1)th
</span><span class="comment">*</span><span class="comment">          columns of VR hold a complex eigenvector, then IFAILR(i) and
</span><span class="comment">*</span><span class="comment">          IFAILR(i+1) are set to the same value.
</span><span class="comment">*</span><span class="comment">          If SIDE = 'L', IFAILR is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          &gt; 0:  if INFO = i, i is the number of eigenvectors which
</span><span class="comment">*</span><span class="comment">                failed to converge; see IFAILL and IFAILR for further
</span><span class="comment">*</span><span class="comment">                details.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Each eigenvector is normalized so that the element of largest
</span><span class="comment">*</span><span class="comment">  magnitude has magnitude 1; here the magnitude of a complex number
</span><span class="comment">*</span><span class="comment">  (x,y) is taken to be |x|+|y|.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            BOTHV, FROMQR, LEFTV, NOINIT, PAIR, RIGHTV
      INTEGER            I, IINFO, K, KL, KLN, KR, KSI, KSR, LDWORK
      REAL               BIGNUM, EPS3, HNORM, SMLNUM, ULP, UNFL, WKI,
     $                   WKR
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.179"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      REAL               <a name="SLAMCH.180"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANHS.180"></a><a href="slanhs.f.html#SLANHS.1">SLANHS</a>
      EXTERNAL           <a name="LSAME.181"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="SLAMCH.181"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANHS.181"></a><a href="slanhs.f.html#SLANHS.1">SLANHS</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="SLAEIN.184"></a><a href="slaein.f.html#SLAEIN.1">SLAEIN</a>, <a name="XERBLA.184"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode and test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      BOTHV = <a name="LSAME.193"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SIDE, <span class="string">'B'</span> )
      RIGHTV = <a name="LSAME.194"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SIDE, <span class="string">'R'</span> ) .OR. BOTHV
      LEFTV = <a name="LSAME.195"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SIDE, <span class="string">'L'</span> ) .OR. BOTHV
<span class="comment">*</span><span class="comment">
</span>      FROMQR = <a name="LSAME.197"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EIGSRC, <span class="string">'Q'</span> )
<span class="comment">*</span><span class="comment">
</span>      NOINIT = <a name="LSAME.199"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( INITV, <span class="string">'N'</span> )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set M to the number of columns required to store the selected
</span><span class="comment">*</span><span class="comment">     eigenvectors, and standardize the array SELECT.
</span><span class="comment">*</span><span class="comment">
</span>      M = 0
      PAIR = .FALSE.
      DO 10 K = 1, N
         IF( PAIR ) THEN
            PAIR = .FALSE.
            SELECT( K ) = .FALSE.
         ELSE
            IF( WI( K ).EQ.ZERO ) THEN
               IF( SELECT( K ) )
     $            M = M + 1
            ELSE
               PAIR = .TRUE.
               IF( SELECT( K ) .OR. SELECT( K+1 ) ) THEN
                  SELECT( K ) = .TRUE.
                  M = M + 2
               END IF
            END IF
         END IF
   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
         INFO = -1
      ELSE IF( .NOT.FROMQR .AND. .NOT.<a name="LSAME.227"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EIGSRC, <span class="string">'N'</span> ) ) THEN
         INFO = -2
      ELSE IF( .NOT.NOINIT .AND. .NOT.<a name="LSAME.229"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( INITV, <span class="string">'U'</span> ) ) THEN
         INFO = -3
      ELSE IF( N.LT.0 ) THEN
         INFO = -5
      ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
         INFO = -11
      ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
         INFO = -13
      ELSE IF( MM.LT.M ) THEN
         INFO = -14
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.243"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SHSEIN.243"></a><a href="shsein.f.html#SHSEIN.1">SHSEIN</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible.
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set machine-dependent constants.
</span><span class="comment">*</span><span class="comment">
</span>      UNFL = <a name="SLAMCH.254"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Safe minimum'</span> )
      ULP = <a name="SLAMCH.255"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Precision'</span> )
      SMLNUM = UNFL*( N / ULP )
      BIGNUM = ( ONE-ULP ) / SMLNUM
<span class="comment">*</span><span class="comment">
</span>      LDWORK = N + 1
<span class="comment">*</span><span class="comment">
</span>      KL = 1
      KLN = 0
      IF( FROMQR ) THEN
         KR = 0
      ELSE
         KR = N
      END IF
      KSR = 1
<span class="comment">*</span><span class="comment">
</span>      DO 120 K = 1, N
         IF( SELECT( K ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute eigenvector(s) corresponding to W(K).
</span><span class="comment">*</span><span class="comment">
</span>            IF( FROMQR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              If affiliation of eigenvalues is known, check whether
</span><span class="comment">*</span><span class="comment">              the matrix splits.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Determine KL and KR such that 1 &lt;= KL &lt;= K &lt;= KR &lt;= N
</span><span class="comment">*</span><span class="comment">              and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
</span><span class="comment">*</span><span class="comment">              KR = N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Then inverse iteration can be performed with the
</span><span class="comment">*</span><span class="comment">              submatrix H(KL:N,KL:N) for a left eigenvector, and with
</span><span class="comment">*</span><span class="comment">              the submatrix H(1:KR,1:KR) for a right eigenvector.
</span><span class="comment">*</span><span class="comment">
</span>               DO 20 I = K, KL + 1, -1
                  IF( H( I, I-1 ).EQ.ZERO )
     $               GO TO 30
   20          CONTINUE
   30          CONTINUE
               KL = I
               IF( K.GT.KR ) THEN
                  DO 40 I = K, N - 1
                     IF( H( I+1, I ).EQ.ZERO )
     $                  GO TO 50
   40             CONTINUE
   50             CONTINUE
                  KR = I
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( KL.NE.KLN ) THEN
               KLN = KL
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
</span><span class="comment">*</span><span class="comment">              has not ben computed before.
</span><span class="comment">*</span><span class="comment">
</span>               HNORM = <a name="SLANHS.310"></a><a href="slanhs.f.html#SLANHS.1">SLANHS</a>( <span class="string">'I'</span>, KR-KL+1, H( KL, KL ), LDH, WORK )
               IF( HNORM.GT.ZERO ) THEN
                  EPS3 = HNORM*ULP
               ELSE
                  EPS3 = SMLNUM
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Perturb eigenvalue if it is close to any previous
</span><span class="comment">*</span><span class="comment">           selected eigenvalues affiliated to the submatrix
</span><span class="comment">*</span><span class="comment">           H(KL:KR,KL:KR). Close roots are modified by EPS3.
</span><span class="comment">*</span><span class="comment">
</span>            WKR = WR( K )
            WKI = WI( K )
   60       CONTINUE
            DO 70 I = K - 1, KL, -1
               IF( SELECT( I ) .AND. ABS( WR( I )-WKR )+
     $             ABS( WI( I )-WKI ).LT.EPS3 ) THEN
                  WKR = WKR + EPS3
                  GO TO 60
               END IF
   70       CONTINUE
            WR( K ) = WKR
<span class="comment">*</span><span class="comment">
</span>            PAIR = WKI.NE.ZERO
            IF( PAIR ) THEN
               KSI = KSR + 1
            ELSE
               KSI = KSR
            END IF
            IF( LEFTV ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute left eigenvector.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLAEIN.344"></a><a href="slaein.f.html#SLAEIN.1">SLAEIN</a>( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
     $                      WKR, WKI, VL( KL, KSR ), VL( KL, KSI ),
     $                      WORK, LDWORK, WORK( N*N+N+1 ), EPS3, SMLNUM,
     $                      BIGNUM, IINFO )
               IF( IINFO.GT.0 ) THEN
                  IF( PAIR ) THEN
                     INFO = INFO + 2
                  ELSE
                     INFO = INFO + 1
                  END IF
                  IFAILL( KSR ) = K
                  IFAILL( KSI ) = K
               ELSE
                  IFAILL( KSR ) = 0
                  IFAILL( KSI ) = 0
               END IF
               DO 80 I = 1, KL - 1
                  VL( I, KSR ) = ZERO
   80          CONTINUE
               IF( PAIR ) THEN
                  DO 90 I = 1, KL - 1
                     VL( I, KSI ) = ZERO
   90             CONTINUE
               END IF
            END IF
            IF( RIGHTV ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute right eigenvector.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLAEIN.373"></a><a href="slaein.f.html#SLAEIN.1">SLAEIN</a>( .TRUE., NOINIT, KR, H, LDH, WKR, WKI,
     $                      VR( 1, KSR ), VR( 1, KSI ), WORK, LDWORK,
     $                      WORK( N*N+N+1 ), EPS3, SMLNUM, BIGNUM,
     $                      IINFO )
               IF( IINFO.GT.0 ) THEN
                  IF( PAIR ) THEN
                     INFO = INFO + 2
                  ELSE
                     INFO = INFO + 1
                  END IF
                  IFAILR( KSR ) = K
                  IFAILR( KSI ) = K
               ELSE
                  IFAILR( KSR ) = 0
                  IFAILR( KSI ) = 0
               END IF
               DO 100 I = KR + 1, N
                  VR( I, KSR ) = ZERO
  100          CONTINUE
               IF( PAIR ) THEN
                  DO 110 I = KR + 1, N
                     VR( I, KSI ) = ZERO
  110             CONTINUE
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( PAIR ) THEN
               KSR = KSR + 2
            ELSE
               KSR = KSR + 1
            END IF
         END IF
  120 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SHSEIN.409"></a><a href="shsein.f.html#SHSEIN.1">SHSEIN</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
